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FINITE ELEMENTS FOR SYMMETRIC TENSORS IN THREE 

DIMENSIONS 

DOUGLAS ARNOLD, GERARD AWANOU, AND RAGNAR WINTHER 



Abstract. We construct finite element subspaces of the space of symme- 
tric tensors with square-intograble divergence on a three-dimensional domain. 
These spaces can be used to approximate the stress field in the classical 
Hellingcr-Rcissncr mixed formulation of the clasticty equations, when stan- 
dard discontinous finite element spaces are used to approximate the displace- 
ment field. These finite element spaces are defined with respect to an arbitrary 
simplicial triangulation of the domain, and there is one for each positive value 
of the polynomial degree used for the displacements. For each degree, these 
provide a stable finite element discretization. The construction of the spaces 
is closely tied to discretizations of the elasticity complex, and can be viewed 
as the three-dimensional analogue of the triangular element family for plane 
elasticity previously proposed by Arnold and Winther. 



1. Introduction 

The classical Lagrange finite element spaces provide natural simplicial finite el- 
ement discretizations of the Sobolev space H^. Similarly, various finite element 
spaces derived in the theory of mixed finite elements, such as the Raviart-Thomas 
and Ncdelec spaces, provide the natural finite element discretizations of the spaces 
if(div) and i?(curl). (These statements are made precise and treated in a uni- 
form framework of the finite element exterior calculus in [3].) In this paper we 
consider the finite element discretization of the space i?(div, f2;§) consisting of 
square-integrable symmetric tensors (or, given a choice of coordinates, symmetric 
matrix fields) with square-integrable divergence. In the classical Hellinger-Reissner 
mixed formulation of the elasticity equations, the stress is sought in If (div, O; S) 
and the displacement in L^(r2;K"). The natural discretization of the latter space 
is evident — piecewise polynomial of some degree without interelement continuity 
constraints — but the development of an appropriate finite element subspace of 
H{div, O; S) to use with these is a long-standing and challenging problem. For 
plane elasticity, the known stable mixed finite element methods have mostly in- 
volved composite elements for the stress [H [TSl [TBI HI]- To avoid these, other 
authors have modified the standard mixed variational formulation of elasticity to a 
formulation that uses general, rather than symmetric, tensors for the stress, with 
the symmetry imposed weakly; see [2 [3 H [HI [H [H [20l [8] . Not until 2002 was 
a stable non-composite finite method for the classical mixed formulation of plane 
elasticity found [10]. This work can be seen as answering the question "what are 



January 16, 2007. 

2000 Mathematics Subject Classification. Primary 65N30, Secondary: 74S05. 



1 



©2007 D. Arnold, G Awanou, R. Winther 



2 



DOUGLAS ARNOLD, GERARD AWANOU, AND RAGNAR WINTHER 



the natural finite element discretizations of i7(div, i7; S)?" in the case of two di- 
mensions. In this paper we address this question for three dimensions. 

Up until recently, there were no mixed finite elements for the Hellinger-Reissner 
formulation in three dimensions known to be stable. In [1], a partial analogue 
of the lowest order element in |10| was proposed and shown to be stable. Here 
we will derive the full analogue of the results of [1^. We construct a family of 
finite element subspaccs of iJ (div, il; §) which, when used to discretize the stress in 
elasticity along with the obvious discontinous piecewise polynomial discretization 
for the displacement, provide stable mixed finite elements for the Hellinger-Reissner 
principle. As in two dimensions, these spaces lead to a finite element subcomplex of 
the elasticity complex, related to it by commuting diagrams. The elasticity complex 
reveals a close connection between the finite clement discretization of i/(div, f2; §) 
and finite element discretization of the space _ff (curl curl*, f2; S), involving a second 
order operator. A key step in our analysis is the identifcation of the conditions 
which arc necessary for a piecewise smooth matrix field to belong to this space. 
These conditions arc derived in Section [5] below. 

We recall the standard mixed formulation for the elasticity equations. Let be 
a contractible polyhedral domain in R^, occupied by a linearly elastic body which is 
clamped on the boundary d^l, and let S and u denote the stress and displacement 
fields engendered by a force / acting on the body. The matrix field S and the vector 
field u can be characterized as the unique critical point of the Hellinger-Reissner 
functional 



over the space H{div,rt;§) x L^(r2;R'^). Here, § is the six dimensional space of 
symmetric matrices and S : T denotes the Frobenius product on §. The given 
compliance tensor A = A{x) : § — > S is symmetric, and bounded and positive 
definite uniformly with respect to x G 51. The divergence operator, div, is applied 
to a matrix field by taking the divergence of each row. Hence, this operator maps 
the space iJ(div,fi;S) into L2(J7;R3)^ 

A mixed finite element method determines an approximate stress field Sh and 
an approximate displacement field as the unique critical point {Sh,Uh) of the 
Hellinger-Reissner functional in a finite element space S/j x V/i C 7J(div, f2; S) x 
i^(51;M'^), where h denotes the mesh size. Equivalently, {Sh,Uh) G S/i x 14 solves 
the saddle point system 



To ensure that the discrete system has a unique solution and that it provides a good 
approximation of the true solution the finite dimensional spaces S/i and Vh must 
satisfy the stability conditions from the theory of mixed finite element methods, 
see [HI [12] . As is well known, see for example [10] , the following two conditions 
are sufficient: 

• divE/i C Vh- 

• There exists a linear operator Hji : H^{fl;E>) Y^h, bounded in C{H^;L^) 
uniformly with respect to h, and such that divH/i S ~ H^divS" for all 
S G iJi(0;S), where : L'^{n;R^) Vh denotes the L^-projection. 

As mentioned above, the construction of finite element spaces which fulfill these 
two conditions has proved to be surprisingly hard. In this paper, we will derive a 
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family of finite element spaces T,h and Vh based on tetrahedral meshes, and show 
that they satisfy these two stability conditions. There is one member of the family 
for each polynomial degree k > 1. The space Vh for the displacements is simply the 
space of all piccewise polynomial vector fields of degree at most k. In the lowest 
order case, fc = 1, the space contains the full space of quadratic polynomials 
on each element, augmented by divergence-free polynomials of degree 3 and 4. 
The local dimension of S/i is 162, or 27 per component of stress on average. The 
analogous space in two dimensions, derived in |10] . was of local dimension 24 (8 per 
component). 

The complexity of the elements may very well limit their practical significance. 
However, we believe that the determination of the natural discretization of space 
_ff (div, Q] S) provides important insight, both into the obstacles to the derivation 
of simpler methods, and for the design of alternative procedures, such as noncon- 
forming methods. 

This paper is organized as follows. After giving some preliminaries remarks 
in Section [2] we present the lowest order element and establish its properties in 
Section [3l A family of higher order elements is then presented in Section |4l Key 
to the analysis of these elements is the description of the polynomial space of 
symmetric matrix fields with vanishing divergence and vanishing normal traces on 
the boundary of a simplex K. The dimension of this space is derived in Section [7| 
based on preliminary results derived in Sections [5] and [6l Furthermore, an explicit 
basis for this space, necessary for the computational procedure, is also given in the 
lowest order cases. Finally, in Section[8]we summarize the results of our construction 
by presenting a discrete analogue of the elasticity complex. 



We begin with some basic notation. If K C M'^ is a tetrahedron, then A2{K) 
denotes the set of the four 2-dimensional faces of K, Ai{K) the set of the six 
1-dimensional edges, and Ao(A') the set of the four vertices. Furthermore, A(A') 
is the set of all subsimplexes of K (of dimension 0, 1, 2 or 3). 

We let M be the space of 3 x 3 real matrices, and S and K the subspaces of 
symmetric and skew symmetric matrices, respectively. The operators sym : M — > § 
and skw : M ^ IK denote the symmetric and skew symmetric parts, respectively. 
Note that an element of the space K can be identified with its axial vector in 
given by the map vec : IK — » : 



i.e., yec~^{v)w = v x w for any vectors v and w. 

For any vector space X, we let X) be the space of square-integrable vector 

fields on fi with values in X. For our purposes, X will usually either be M, M.^, 
or M, or some subspace of one of these. In the case A = R, we will simply write 
L^(ri). The corresponding Sobolev space of order k, i.e., the subspace of L'^{V,]X) 
consisting of functions with all partial derivatives of order less than or equal to k 
in L'^{n;X), is denoted and its norm by || • ||fc. The space i7(div,rj;S) 

is defined by 



2. Notation and preliminaries 




H{div,n;S) = {Te L^{n;S) \ divT G L^{n;R^)}, 
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where the divergence of a matrix field is ithe vector field obtained by applying the 
divergence operator row-wise. For a vector field w : f2 — > R'^, gradw is the matrix 
field with rows the gradient of each component, and the symmetric gradient, £(v), 
is given by e{v) = symgrad?;. Furthermore, 

(93^2 - 92^3 \ 
—d^vi + 9if3 
92^1 - dlV2 / 

If we consider a linear coordinate transformation of the form x = Bx + h, with 
corresponding vector fields v and v related by v{x) = {B')~^v(x), then we have 

skwgrad^u = (i?')^^(skwgrad^ u)_B~^. 

Here B' denotes the transpose of B. In particular, if B is orthogonal, i.e., B'B = /, 
then V = Bv and 

(2.1) skwgrad^ v = i3(skwgrad^ v)B' . 

As for the divergence and the gradient operator, the operator curl acts on a 
matrix field by applying the ordinary curl operator to each row of the matrix. 
The operator curl* is the corresponding operator obtained by taking the curl of 
each column. Alternatively, we have curl* T = (curlT')' for any matrix field T. 
The second order operator curl curl* maps symmetric matrix fields into symmetric 
matrix fields. Let S : M ^ M be the algebraic operator ET = T' - tr(T)/, where 
/ is the identity matrix. Then S is invertible with T = T - tr(r)//2. The 
following identities are useful: 

(2.2) vec skw curl T = -i div ST, reC°°(n,M), 

(2.3) curlT = SgradvecT, TeC°°(r2,K), 

(2.4) trcurlT = -2divvccskwT, reC°°(17,M). 

These formulas can be verified directly, but they are also consequences of the dis- 
cussions given in [SI [3], cf. Section 4 of [S] or Section 11 of 

For if C M'^ we let Vk{K;X) be the space of polynomials of degree fc, defined 
on K and with values in X. We write Vk or VkiK) for 'Pk{K;R). The de Rliam 
complex has a polynomial analogue of the form 

(2.5) R ^ Vk+3 ^ Vk+2{K;R') ^ rk+,{K;R') ^ Vk ^ 0. 

In fact, this complex is an exact sequence [5]. 

In recent years differential complexes have come to play a significant role in the 
design of mixed finite element methods [31 [TOl [51 [5] . For the equations of elasticity, 
the relevant differential complex is the elasticity complex. In three space dimensions, 
the elasticity complex takes the form 

where T is the six-dimensional space of infinitesimal rigid motions, i.e., the space 
of linear polynomial functions of the form x i-^ a + b x x ior some a, 6 G M'^. It is 
straightforward to verify that the elasticity complex is a complex, i.e., the compo- 
sition of two successive operators is zero. In fact, if the domain CI is contractible, 
then the elasticity complex is an exact sequence; see [SI [9]. 
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An analogous complex with less smoothness is 
T H\n; R3) ^ //(curlcurl* , 17; S) (div, il; S) ^ L^i^; R^) _ q, 

where Jf(curlcurr , 17; §) = {5 £ L^{n;S) \ curl curl* 5 e L^{n-S)}. 

There is also a polynomial analogue of the elasticity complex. Let K C M'' be 
tetrahedron and fc > 0. The polynomial elasticity complex is given by 

(2.6) T ^ Vk+4{K;M.^) ^ Pk+3{K;S) Vk+i{K;S) ^ rk{K;R^) ^ 0. 

This complex is an exact sequence. To prove the exactness, we first show that if 
5 is a matrix field in Vk+si^]^), with curl curl* 5' = 0, then S = e(u) for u = 
(ui,U2,U3)'^ G 7'fc+4(^^;R^)• Clearly S = e{u) for u e C°°{K;'&.^). It is enough to 
show that all second derivatives of u^, fc = 1, 2, 3, are in ■Pfe+2(S7;M). This follows 
from the identity dijUk = die{u)jk + dj e{u)ik — dke{u)ij. Next, we show that 
if S* e 7'fc+i(17;§), and divS" = 0, then S = curl curl* T for some T e 7'fc+3(r2;S). 
First we observe that since div S" = it follows from the fact that (|2.5p is exact that 
S = curl?7 for some U € 'Pk+2{K; M). Furthermore, since S is symmetric it follows 
from (|2.2p that div S ?7 = 0, and as a consequence, using (|2.5p once more, we obtain 
that EU = curlT for some T G Vk+3{K;M), or S = curlJJ = curl S^^ curl T. 
However, by (|2.3[) we have curl curl skwT = curl grad vec skw T = 0. Hence, 
we can take T G Vk+siK^S). Finally, we observe that if T is symmetric, then 
(|2.4p implies that trcurlT = 0, and therefore S = curl curl T = curl curl* T. 
To establish the surjectivity of the last map, one can use the fact that dim Vk = 
(fc + l)(fc + 2)(fc + 3)/6 to verify that the alternating sum of the dimensions of the 
spaces in the sequence is zero. The same arguments show that p.6p is exact for 
k = —1, —2, or —3, if Vj is interpreted as the zero space for j < 0. 

Let denote a family of triangulations of by tetrahedra with diameter 

bounded by h. We assume that the intersection of any two tetrahedra in 7^ is 
either empty or a common subsimplex of each. The family {7^} is also assumed 
to be shape regular in the sense that the ratio of the radii of the circumscribed 
and inscribed spheres of all the tetrahedra can be bounded by a fixed constant. 
Furthermore, we will use the notation Aj{Th), for j = 0, 1, 2, to denote the set of 
vertices, edges, and faces, respectively, associated with the mesh Th. In Section|3]we 
will define a family of finite element spaces 'Eh C iJ(div, 12; S) and Vh C i^(17; R) for 
the elasticity problem consisting of piecewise polynomial spaces with repect to Th of 
arbitrarily high polynomial order. However, we will first consider the lowest order 
case of this family in Section [3] below. All our spaces will have the property that 
divS/i C Vh- Furthermore, we will identify a corresponding projection operator 
H/i : H^{fl;§) E/i satisfying the commutativity relation 

(2.7) divH/iT = H^divT, T e H\il;S), 
and the bound 

(2.8) ||n,T||o<C||T||i, TeH\n;§), 

with constant C independent of h. Here H^ : L^(r2;K'^) Vh is the projec- 
tion. It is a consequence of the general error bounds derived in [M], cf. also [lO] . 
that the properties above imply that (E/i, 14) is a stable pair of elements for the 
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discretization (jl.ip . and that the error bounds 

(2.9) ||5-5^||o < ||(/-n^)5||o 

(2.10) \\u - u^llo < ||(/ - TlXMo + 4(1 - n^)5||o 

holds, with a constant c independent of h. Here {S,u) is the unique critical point 
of the Hellinger-Reissner functional over iJ(div, f2;§) x L^(f2;M^) and {Sh,Uh) G 
S/i X Vh the corresponding finite element solution. In addition, div5?i = Il^divS'. 

3. The lowest order element 

We first describe the restriction of the lowest order spaces "Eh and Vh to a single 
tetrahedron K ^ T^. Define 

Ex = {T G Vi{K;S) I divT G PiiK^W") }, Vk = Vi{K;R^). 

The space Vk has dimension 12 and a complete set of degrees of freedom is given by 
the zero and first order moments with respect to K . The space has dimension at 
least 162 since the dimension of Vi{K\ S) is 210 and the condition div T G Vi {K; R^) 
represents 48 linear constraints. We will show that diniE^ = 162 by exhibiting 
162 degrees of freedom which determine the elements uniquely. Define 

M{K)^{T e ViiK; §) | div T = 0, Tn = on dK }. 

It will follow from Theorem 17.21 below that the dimension of Ai{K) is 6. Further- 
more, in Section [7] we will also give an explicit basis for this space. Using this basis, 
we can state the 162 degrees of freedom for the space S^. 

Lemma 3.1. A matrix fieldT G Yik is uniquely determined by the following degrees 
of freedom 

(1) the values of T at the vertices of K , 4 x 6 = 24 degrees of freedom, 

(2) for each edge e G /S.i{K) with unit tangent vector s and linearly independent 
normal vectors n_ and the constant, linear and quadratic moments over 
e of s'Tn^, s'Tn^, n'_Tn^, n'_^Tn^, n'_Tn-^-, 6 x 3 x 5 = 90 degrees of 
freedom, 

(3) for each face f G A2{K), with normal n, the constant and linear moments 
over f ofTn, 4 x 3 x 3 = 36 degrees of freedom, 

(4) the average of T over K , 6 degrees of freedom, 

(5) the value of the moments J^^T -.U dx, U G A4{K), 6 degrees of freedom. 

Proof. We assume that all degrees of freedom vanish and show that T = 0. Since 
T = at the vertices, the second set of degrees of freedom imply that Tn = on 
each edge for both faces meeting the edge. By the third set of degrees of freedom 
we obtain that Tn = on each face of K. For v = divT G Vi{K,M.^) we have 

/ dx ~ — T : e(u) dx + Tn ■ v dxj = — T : ev dx ~ 
J K JK JdK Jk 

by the fourth set of degrees of freedom. Here and below, dxf denotes the surface 

measure on dK. We conclude that divT = 0, and, by the last set of degrees of 

freedom, that T 0. □ 

We now describe the finite element spaces on the triangulation T^. We denote 
by Vh the space of vector fields which belong to Vk for each K G Th and by T,h 
the space of matrix fields which belong piecewise to T,k, and with the continuity 
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conditions induced by the degrees of freedom. In particular, for T e S^, the 
normal components Tn are continuous across all faces / G A2(7/i) and, hence, 
S/i C i?(div, 51; §). In addition, if T G E/i, e G Ai(7^), and s and n arc vectors 
which are tangential and normal to e, respectively, then s'Tn is continuous on e. 

It remains to define an interpolation operator II/j : H^{fl; S) S;, which satisfies 
p.7p and (|2.8p . Because of the vertex and edge degrees of freedom, the canonical 
interpolation operator for E/i, Hj^, defined directly from the degrees of freedom, 
is not bounded on i7^(ri;§). In order to overcome this difficulty we introduce the 
operator 11° : H^{fl; §) — > E/j defined from the degrees of freedom above, but where 
the vertex and edge degrees of freedom are set equal to zero, i.e., we have 

X e Ao(T„), 

e e Ai{Th), V £ ■P2(e;M^), n e e^, 

K e n, 

K G Th, U G MiK). 
The commutativity property (|2.7p for IlJJ follows from ()3.3p and (|3.4p since 
/ (divn$Jr-divT)-wdT = - / {UlT-T) : e{v)dx+ {UlT~T)n-vdxf = 0. 

Jk JK JdK 

The uniform boundedness (j2.8[) can be seen from a standard scaling argument 
using the matrix Piola transform. Let if be a fixed reference tetrahedron and 
F = Fk '■ K K be an afRne isomorphism of the form Fx = Bx + b. Given 
a matrix field T : K ^ §>, define T : K ^ S hy the matrix Piola transform 
T{x) = Bf{x)B^, with x = Fx. Using divr(x) = Bdivf (i), it is easy to verify 
that T G if and only if T G E^. Furthermore, as in [11 [10] a scaling argument 
can be used to verify the uniform boundedness condition (|2.8p for the operator 
n°. We can therefore conclude that the operator 11° satisfies the two conditions 
p.yp and (|2.8p . However, the operator 11° lacks good approximation properties. 
Therefore, in order to obtain error estimates from the general bounds (|2.9p and 
(|2.10p a more accurate interpolation operator is needed. 

Consider the modified interpolation operator Hh : H^{il; §) — > 'Sh of the form 

(3.6) nh=nl{I -Rh) + Rh, 

where Rh ■ L'^{fl;§) E/i is the Clement operator onto the continuous piccewise 
quadratic subspace of E^ |13| . This operator satifics the bounds 

\\RhT-T\\,<ch"'-mT\U 0<j<l, j<m<3. 

As a consequence of this bound, and the boundedness (|2.8p of 11°, we obtain the 
estimate 

(3.7) ||n;,r-T||o<c/i"||r||™, l<m<3 

for the interpolation error. Furthermore, since satisfies I^Ji) and div Rh = 
divi?ft, we conclude that Ilh satisfies p.7p . 



(3.1) 


n°r(x) 


= 0, 


(3.2) 


J nlTn ■ V ds 


= 0, 


(3.3) 


j(T-nlT)n-vdxf 


= 0, 


(3.4) 


f {T-UlT)dx 


= 0, 


(3.5) 


f {T-UlT):Udx 


-0, 
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We also recall that the projection operator : L'^{Q;R^) Vh satisfies the 
error estimate 

(3.8) \\UXv^v\\o<ch"'\\v\U, 0<m<2. 

The estimates p.7p and p.Sp . combined with the basic error bounds (j2.9p and 
p.lOp . and the fact that div Sh = ti^ div S, imply the following error estimates for 
the finite clement method generated by E/j x V^,. 

Theorem 3.2. Let (S,u) denote the unique critical point of the Hellinger-Reissner 
functional over i7(div, $7; §) x £^(51; M?) and let (Sh, u^) be the unique critical point 
over S/j X Vh ■ Then 

\\S — Sh\\a 
II div 5 - divS';i||o 
ll^t - Uh\\o 

Remark. In general it seems not possible to lower the polynomial degree of the 
stress space E/i introduced above. However, as in [10| . a minor simplification is 
possible. On each tetrahedron K Th we take the restricted displacement space 
Vk to be the rigid motions T C 7-'i(/^;R'^) and the corresponding stress space to 
be 

Ea' = {T G V4{K;S) I divT e T}. 

Clearly dim E^- > 210 — (60 — 6) = 156. In fact, dimE^ = 156 and a complete set of 
degrees of freedom is obtained by removing the six average values of T represented 
by (4) in Lemma 13.11 The proof of the fact that these degrees of freedom are 
unisolvcnt for E^" follows by a simple modification of the proof of Lemma 13 . 1 1 above . 
Just observe that if v = divT G T then e(u) = 0. However, the simplified clement 
is less accurate, since the stress space lacks some quadratics, and the displacement 
space some linears. Instead of the error estimates given in Theorem 13.21 we obtain 
at most 0{h'^) convergence for US' — 5/j||o, and at most first order convergence for 
I div(5 - 5*^)110 and \ \u - UhWo- 

4. A Family of Higher Order Elements 

In this section we describe a family of stable element pairs, one for each degree 
fc > 1. The lowest order case fc = 1 is the one treated above. We first describe the 
elements on a single tetrahedron. Define 

Ek = { T e Vk+^iK; §) I divT G Vk{K- M^) }, Vk = Vk{K; W"). 

Then 

k + 3\ _ (fc + 3)(fc + 2)(fc + 1) 
3 J ~ 2 ' 

dimE^ >dk: = dimVk+siK;^) - [dimVk+2{T;M.^) - dimVkiTiR^)] 

= 6(^;*^)-3('-;^)+3(^;^)=fc^ + m^ + 56fc + 93. 

Notice that the space e[7'fc(A', M^)] has dimension (fc + 3)(fc + 2)(fc + l)/2 - 6. 
Analoguous to the lowest order case, we define the space 

Mk{K) = { T G Pk{K; §) I divT = 0, Tn = on dK }. 



< c/i"||S'|U, 1<TO<3, 

< c/i" II div 5||,„, 0<m<2, 

< c/i™||u|L+i, 1 < m < 2. 



dim Vk = 3 
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We will prove in Section[7l TheoremESl that dimMkiK) is (fc + 2)(fc-2)(fc-3)/2 
for fc > 4. 

The degrees of freedom for Vk are the moments of degree less than or equal to 
k with respect to K. A unisolvcnt set of degrees of freedom for Y.k arc given by 

(1) the values of T at the vertices of isT, 4 x 6 = 24 degrees of freedom, 

(2) for each edge e £ Ai{K) with unit tangent vector s and linearly indepen- 
dent normal vectors n- and n+, the moments of degree at most k + I over 
e of s'Tn_, s'Tn+, n'_Tn., n'^Tn+, n'_Tn+, 6 x (fc + 2) x 5 = 30fc + 60 
degrees of freedom, 

(3) for each / S A2{K) with normal vector n, the moments of degree at most fc 
over / for Tn, 3 X 4 X (fc + 2) (fc + l)/2 = 6fc2 + 18fc + 12 degrees of freedom, 

(4) T : C/rfx, [/ e eiVx), (fc + 3)(fc + 2)(fc + l)/2 - 6 degrees of freedom, 

(5) Jj^T -.Udx, U e Mk+3{K), (fc + 5)(fc + l)fc/2 degrees of freedom. 

The proof that this set of functionals is unisolvent for the space is almost 
identical to the lowest order case, and it is easily checked that their numbers sum 
up to dk- Hence, we have shown that dimE^ ~ dk- Furthermore, in Section [7] we 
will give an explicit basis for the space Mk{K) when fc = 4 and fc = 5. 

The finite element space Vh C i^(f^; R^) consists of all vector fields which belong 
to VkiK; R'^) for each K G T/^, while the corresponding stress space E/j is the space 
of matrix fields which belong piecewise to T,f^, and with the continuity conditions 
induced by the degrees of freedom. In particular, this implies that the normal 
components Tn, for T G E/i, are continuous over all faces in A2{Th). Hence, as in 
the lowest order case we have that E/i C H{div, fl; S). 

The projection H^ onto Vh satisfies the estimate 

(4.1) \\IlXv~v\\o<ch"'\\vln, 0<m<fc + l. 

We also introduce the Clement interpolant Rh : L^(ri;§) E/i defined as the §- 
valued version of the standard scalar Clement interpolant into continuous piecewise 
polynomials of order fc + 1. Hence, the operator satisfies 

\\RhT - T\\, < ch"'~^\\T\U < J < 1, J < m < fc + 2. 

Furthermore, we define the modified canonical interpolation operator H/j by (|3.6p . 
where the operator H° is defined in complete analogy with the lowest order case, 
by setting the degrees freedom associated with Ao(7^) and Ai{Th) equal to zero. 
Then the operator H/i satisfies (|2.7p and (|2.8p . and the error bound 

(4.2) ||H,,T-rl|o<c/i"||r|j„„ l<m<fc + 2. 

As above, the interpolation estimates (|4.ip and (|4.2|) . and the error bounds p.9|) 
and (j2.10p . leads to the following error estimates. 

Theorem 4.1. Let (S,u) denote the unique critical point of the Helling er-Reissner 
funtional over H{div,n;S) x L^in^m^) and let (Sh.uh) he the unique critical point 
over E^ x V/j. Then 

llS'-S'ftllo < c/i™||5||,„, 1 < m < fc + 2, 
||divS'-divS',,||o < c/i^lldivS-ll™, < m < fc + 1, 
\\u - Uhh < c/i"||u||„+i, 1 < m < fc + 1. 
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5. Some properties of vector fields and matrix fields 

In order to perform a more detailed analysis of the spaces A4k{K) and the finite 
element spaces S^, we will need some basic properties of vector fields and matrix 
fields. These properties will be reviewed in the present section. 

5.1. Identities for vector fields. Let n denote a fixed unit vector in R^, P„ = 
nn' the orthogonal projection onto Mn, / = the plane orthogonal to n, and 
Qn = I — Pn is the orthogonal projection onto /. Furthermore, set 

(0 
-ris ni 
712 -ni 

so that CnV = V X n. The following identities are easily checked: 

For any vector field v = (wi, V2, v^)' in R^, we obviously have v = PnV + QnV and 
curlw ~ cuiiP„7;+curl(5„w ~ P„ curl P„7;+P„ curl QnW+Qn curl P„7>+(5„ cwclQnV. 
It is elementary to verify that 

dv 

Pn curl PnV ^0, Qn curl QnV = C„ — , 

077 

if 71 = 63 = (0, 0, 1)'. In view of the transformation formula (|2.ip . these identities 
hold for an arbitrary unit vector 77. Furthermore, we define 

(5.1) TotfV = P„ curl 71 = P„curl(5„7J ~ — (div C„7;)7i, cwclfV = Q„curlP„7;. 
With this notation we obtain the decomposition 

dv 

(5.2) curl 71 = rot f 7) + curl f 7) + C„ — — . 

dn 

We also define the tangential gradient grady </> = Qn grad0, for a scalar field (f>. For 
a vector field v, grad^T; = (grad v)Qn is the matrix field with rows equal to the 
tangential gradients of the components of v, and we let 

*^/(^) = ^{grad/((3„7;) + [grad/((3„7))]'} = Qn £(«)<?« 

be the tangential part of the symmetric gradient. Note that the definitions of 
gradjT;, £/(?;), curl/?;, and rot/7; do not depend on the choice of the unit normal 
rt to /. The identity 

(5.3) curl/7) = — C„ grad(77'7;) = — C„ grad/(77'7;) 

can be easily verified in the special case 77 = 63, and holds in general. 

5.2. Identities for matrix fields. We extend the operators curl, curl/ and rot/ 
to act on and yield 3x3 matrix fields by applying the vector operations row-wise. 
More precisely, rot/ 5 = (curlS')P„ = (curl5'Q„)P„ and curl/ 5 = (curl 5Pr,)(5n- 
We notice that, for any constant matrix A, curl AS* = A curl S*. We also recall 
that curbs' = (curl 5")' is the corresponding operator obtained by applying the 
curl operation to each column. The corresponding column operators rot/ and 
curl/ are defined similarly, i.e., rot/ 5* = P„ curl* 5 = P„ curl* QnS and curl/ S = 
Qn curl* PnS. For a given row vector 7;, grad/ v — (grad/ v')' = Qnig^^iv')' , which 
is the matrix whose columns are the tangential gradients of the components of 7;. 
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We now extend the decomposition (|5.2p to curl and curl* . It is easy to see that 
CnS results in C„ applied to the columns of S, while — 5C„ is C„ applied row-wise. 
It follows that 

dS 

(5.4) curls' = curl/ 5 + rot/ S' - —C„, 



and 



dn 

dS 



(5.5) curl* S = curl*. S' + rot*f S + C„ ^ , 

•' on 

where dS/dn is obtained by taking the directional derivative of each component. 
Furthermore, the identities 

(5.6) curl/ S* = (grad/ 5ri)C,i, curl/ 5 = —C,i grad/ n'S* 

are just matrix analogues of (|5.3p . Note also that from the definitions of the 
operators rot / and rot ^ we have 

(5.7) P„(curlcurr S)Pn = P„(rot/ curl* S) = rot/(P„ curl* S) 

(5.8) = rot/ rot} 5 = rot/ rot}(Q„S'(5„). 

We will also need exact sequences relating spaces of functions defined on a two 
dimensional space. Let / = . In analogy with (j2.6p the following two dimensional 
complexes are exact: 

n e f rot f rot f 

(5.9) T/ Pk+3i.f; 0„m3) M Vk+2{f; QnSQn) Vk{f;RPn) ^ 0, 

(5.10) 

Vl{f;R) ^ rk+3{f;M.) ^^^^^^^ Vk+liflQu^Qn) ^ Pk{f;QnSPn) ^ 0. 

Here, T/ is the 3-dimensional space of vector fields on / of the form v{x) ~ Qnw{x) 
for some w £T. 

5.3. Integration by parts. Above we discussed the operators P„, Qm and C„ 
with respect to a fixed linear space / with a unit normal vector n. If is a bounded 
subset of R^, with a piecewise smooth boundary Qfi, we define these operators on 
5fi with respect to the tangent space / = f[x) and the outer unit normal vector 
n = n{x) for x £ dVL. Furthermore, the tangential differential operators grad/, 
curl/, rot/, e/ and so on are then defined pointwise on dO.. 
If w is a smooth vector field on f2, then 

/ V ■ curl If dx ^ curlw ■ w dx ~ I C„u ■ w dxf, 
Jn Jn Jan 

where n denotes the outward unit normal to O, and dxf is the surface measure on 

5r2. Moreover, if v and w are smooth vector fields on dU, and / is a piecewise 

smooth submanifold, then we have 



rotfv-w dxf ~ I v-cm\fW dxf — / {s'v)(n'w) ds, 
'f Jf Jdf 

where ds denotes the arc length measure on df. 

We next extend the previous integration by parts formula to matrix fields. Let 
S and T be two smooth matrix fields on M"^ not necessarily symmetric. We have 

(5.11) / : curlT / curlS':rda;+ / SCn'-T dxf 
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and 

(5.12) [ rot/ S:T dxf= [ S : cuAjT dxf - [ {Ss) ■ (Tn) ds. 
Jf Jf Jdf 

Given a symmetric matrix field 5* define A/(S') : / Qn^Qn by 

wfiere dnS := dS/dn. Hence, Aj^(iS') is a symmetric matrix field defined on /. If 
T = e{v), where w is a vector field then we have 

2ef{Tn) = grady grad^(?7/i;) + Qndn<^{v)Qn- 

Hence, we obtain that 

(5.13) A/(e(i;)) = grader grady(n'i;). 

The tangential-normal components of the matrix field curl curl* S on f can be 
expressed in terms of Af{S). Indeed, by the definition of the operator rot/ and 



()5.5p we have 

C„ (curl curl* S)Pn = C„ rot/ curl* S' = rot/ C„(curl/ S + rot/ S + CndnS). 

However, d rot/ S* = and, by (|5.5p . Cn curl/ 5* = Qn grad/ n'S. Hence, 

(5.14) C„(curlcurl* S)Pn = rot/ Q„(grad} n'S - dnS)Qn = rot/ A/(S'). 

The next lemma indicates how the operator A/ arises when integrating curl curl* 
by parts. 

Lemma 5.1. Let S and T be two smooth matrix fields, with S symmetric, on Q, 
where fl is a bounded subset o/R'^ with a piecewise smooth boundary 9f2. Then 

S : curl curl* T dx ~ ( curl* curlS" : T dx — ( CnSCn : dxf 



[ [Cn rot/ 5 - (rot} S)Cn + C„A/(5)C„] -.Tdxf. 
Jdn 



I an 

Proof. Since Sfi has no boundary, it follows from (|5.12p that 



(5.15) 



/ rot f S : T dx f ^ / S : curl / T dxf. 
Jan Jan 



Then using (|5.4p and the antisymmetry of C„ wc obtain 

f S : curl T dxf = f S : curl fTdxf+ f S : rot fTdxf—f S : — — C„ dx / 
Jan Jan Jan Jan on 

(5.16) =/ (rotf +cnTlf)S :Tdxf + f SCn-^dxf. 

(So far we have not used the assumption that S is symmetric.) Next, from (|5.1ip 
we get 

/ S ■ curl curl* T da; = / curl 5 : (curir')'da; + / SCn : (curlT')'da; 
Jn Jn Jan 

= [ (curl 5)' : curlT'dx - [ C„S" : curlT'dx/ 
Jn Jan 
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and 



/ (curls')' : curl T'dx [ cml{cur\ S)' : T' dx + [ (curl S')'C„ : T' dx/. 
Jn Jn Jan 



In Jn Jan 

Since [curl(curl S*)']' = curl* curlS" we obtain the identity 

curl* curl S : T dx 



n 



I : curl curl* T da- + / C„ curl S" : T da;/ + / C„S' : curl T' da;/. 
Jn Jan Jan 



In Jan Jan 

Using (|5.16|) . with S replaced by CnS, we can rewrite the last term as 

f if if dT' 

I CnS : curlT' dx f = / C„(curl/ + rot/)S' : T' da;/ + / CnSCn ■ — dxf 
Jan Jan Jan c'" 

= - / [(curl} + rot}) S']C„ : Tdxf + f C„5C„ : ^ da;/. 

Thus we obtain 

(5.17) / S : curl curl* T dx ^ [ curl* curlS" : T dx — [ CnSCn ■ -rr- dxf 
Jn Jn Jan dn 

-[ CnCmlS -.Tdxf + [ [(curl}+rot})5]C„ : Tda;/. 
Jan Jan 

In order to see that (j5.17p is equivalent to the desired identity it is enough to show 
that 

(5.18) C„ curl S - (curl} S')C„ = C„A/ (S')C„ + C„ rot/ S. 

However, from (|5.6p it follows that C„ curl/ S — (curl} S)Cn — 2Cn e{Sn)Cn, and 
hence ([5718]) follows from ([O]) . □ 

The lemma can be used to determine when a symmetric matrix field S 011 f2 
which is piecewise smooth with respect to a given triangulation belongs to the 
space 7J(curlcurr, ri; §). This holds if and only if curl curl* 5 defined piecewise 
coincides with curl curl* S defined in the sense of distributions. Integrating against 
a smooth test function and using the lemma, we find that a necessary condition is 
that CnSCn, or equivalently QnSQn, is continuous across element faces. In this 
case the quantities C„rot/5 and (rot} S')C„ are also continuous. Again invoking 
the lemma again, we obtain the following result. 

Theorem 5.2. Assume that S £ L^(ri;S) is piecewise smooth with respect to a 
triangulation T of Q. Then S G i7(curlcurr, fi; S) if and only if QnSQn o.nd 
J^f{S) are continuous across each face in A2(T). 

6. Polynomial matrix fields on a single tetrahedron 

In order to study the space S/i C iJ(div, $7; §) introduced above, we will need 
to study finite element subspaces 8/1 of iJ(curlcurr, 17; §). The discussion in the 
present section will be restricted to a single tetrahedron K , but based on the results 
derived here we will define the space Qt in the final section of the paper. In fact, 
we will present a complete discrete elasticity complex of the form 

TTi Mir r\ curl curl* „ div _ 
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where Wh C i?^(fi;R'^) and Qt C (curl curl*, 17; §) are piecewise polynomial 
spaces with respect to the triangulation Th- 

Let K C M'' be a fixed tetrahedron and define the polynomial space 

(6.1) = Afk{K) = {S€ VkiK; §) | Q„5g„ = Af{S) = 0, / G A2{K)}. 

Most of the discussion in this section is devoted to computing the dimension and a 
basis for this space. However, first we need some additional notation. 

If / G A2{K), we denote by ft,/ the perpendicular distance from the opposite 
vertex to / and by n = Uf the outward normal vector to /. If e is an edge, we 
let s — Se denote one of the unit vectors parallel to e. When the edge e belongs 
to the face /, we write m = m^j for the unit vector in /, normal to e, pointing 
from e into /. See Figured] When the notations /+ and /_ arc used to denote two 
faces, the corresponding normals will be denoted n+ and 7i_, and the perpendicular 
distances and respectively. The notations m_|_ and m_ will also be used to 
denote rrie.f^ and mej_ where e is the edge common to /_)_ and /_. 




Figure 1. The [uf , Sf-^m^.f) coordinate system for a face / and 
edge e of the tetrahedron K. 

The barycentric coordinates on K will be labelled by the faces. That is, they 
are \f £ Vi{K\K) determined by Ay = on / and X)/-^/ = 1 on K. We recall 
that grad A/ = —Uf/hj. Let g G A.{K) be a face of dimension m with vertices 
Xi„,Xi^, . . . , . For TO = 0, g is a vertex, for to = 1, g is an edge, and so on. We 
define the bubble functions bg = \f.^\f._^ ■■■Xf.^^, where fi^, is the face opposite 

vertex . For d > 0, (K) = span{ A]:°^ A;^'^ ' " '^^f" I -^o H h jm = ^ }, so that 

dimV^iK) = dimPd{g). Note that if x is a vertex opposite face /, then bx = A/ 
and PS{K) = RXj, while bx = 11/ -^Z and = Vd- For a given face /o G A2{K), 
bfo = bx/Xfo = Uf^fo -^/' and for a given edge e G Ai{K), be = bK/{Xf_Xf^), 
where /_ and /_|- are the faces containing e. 

The monomials of degree k in the barycentric coordinates form a basis for VkiK), 
and by grouping together terms according to which coordinates enter the monomial, 
we can uniquely represent any p G "Pfe as 

(6-2) E P.en'-i-dim,W- 

geA{K) 

The standard Lagrangian degrees of freedom for p E Pk (K) are the values of p at 
the vertices, the moments of p of degree at most — 2 on each of the edges of K, the 
moments of degree at most fc — 3 on each of the faces, and the moments of degree 
at most k^A on K . From the vertex values of p we may determine the polynomials 
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Pg in (|6.2p for g G Ao(/^). From these and the edge moments we may determine as 
well the Pg for g e Ai{K), etc. 

Of course analogous considerations apply to Vk{K;X) for X a vector space. 
In particular, we have the representation p — J2g£A{K)^aPa ^^'^ P ^ T-'kiK-jX) 
where now pg G ^fc_i_dimg(-^' ■^)^ ^^^'^ space of X- valued polynomials on K whose 
components with respect to a basis of X belong to T^k-i-dimgi-^)- 

If e = /_ n /+ e Ai(A'), with f-,f+e A2iK), we let Ge e S be the matrix 

We note that GgS = and m'_Gem^ ~ m'^Gem^ ~ 0. 
Lemma 6.1. For k > 0, the dimension of the space 

Af^ = J^^iK) :^{Se Vk{K- S) | Q„5g„ = 0, / e ^^{K) } 
is (k + l)k{k ~ 1). 

Proof. For a subsimplcx g G A(/v), let us first define 

N'^ = { S* G § I QnSQn = for all faces containing g }. 

Clearly, ii g = K then dim = 6, if g G A2{K) then dim A^f = 3, and if 
g G Ao{K) then dim A^^ = 0. Finally, if 5 G Ai{K), i.e., 17 = e is an edge, then 
dim = 1. In fact, the space N"^ is then spanned by the matrix Ge introduced 
above. 

If 5 G J\fj! then from (|6.2|) we obtain the representation 
(6-3) 5= E ^oSs, SgeV',_,_^,^^{K;N^). 

As a consequence 
dim<= ^ dim7't,_,,„^(A-7V^) 

= 6{k - 1) + 6(fc - l)(fc - 2) + (fc - l)(fc - 2)(fc - 3) = (fc + l)fc(fc - 1). □ 
Before we arc able to compute diniA/fc, we need to establish several lemmas. 
Lemma 6.2. If S £ Afk, S is zero on each edge. 

Proof. As above let e = /_ n /+ G Ai(A'), with /+ G A2{K). Then 

m'_n^ = m'_^_n^ < 0. 

In fact, the transformation 

is a rotation in the plane orthogonal to e. 

Let S G A/l-- Since A/l- C A/]? we know that S\e = pGe, where p G Vkie). 
Therefore, if we can show that 

(6.5) m'^Sn^ + m'_Sn^ = 

on e, then p{m'_^_n^ + m'_n-^) = 2pm'_^n^ = 0, and, as a consequence, 5* is zero on 
e. It therefore suffices to show (|6.5p . 

On e, we must have s'Af{S)m = 0, i.e., 

(9s(m'57i) + d^is'Sn) - dn{s' Sm) = 0, 
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where / is either /_ or /+. By adding this property for the two faces we obtain 

— ds{m'^Sn+ + m'_Sn-) 

= [drags' Sn^) ~ 9„_(s'5m_)] - [c'„Js'5m+) - 

However, the right hand side here is zero as a consequence of the fact that the 
transformation (|6.4[) is a rotation. In fact, this property imphes that 

9„_(i;'7i_) - dn_{v'm^) = dn+{v'm+) - dm^{v'n+) 

for any smooth vector field v on K. Hence, we can conclude that m^S'n-|_ + m'_S'n_ 
is a constant along e. and since it is zero at the vertices, (|6.5p holds. □ 

Lemma 6.3. Let 

/6A2(K) 

Af(U)\e — for each face f and each edge e of f}. 
Then diuiAfk.dK = 6(P - 6fc + 10). 

Proof IW ==J2fbf Uf , then U G JV° if and only if each coefHcient Uf eVi_^{K; S) 
satisfies 

QnUfQn = on/. 

Hence, this property is assumed to hold. We have A/(C/) = on an edge e C / if 
and only if the three terms s'Af{U)s, s'Af(U)m and m' Kf{U)m vanish there. For 
any fixed unit vector t and e G /S.i[K), we have 

(6.6) dtU = Y,{dtbs)Us = -heC^Uf„ + ^^Uu) one, 

where /_ and /+ are the two faces meeting the edge e and we have used that 
grad6/_ = —n+be/h+ on e. 

Recall that s'Af{U)s = 2d,{s'Un) - 9„(s'C/s). Since L/ = on e, 

s'Af{U)s = -dnis'Us) on e. 

However, since QnUfQn = on the face /, we have s'Uf_s = s'Uf^s = on e. By 
(|6.6p . with t = n, we conclude that s'Af{U)s = on e. 

Next, similar considerations for s'Af{U)m = ds{m'Un) + dm{s'Un) — dn{s'Um) 
give 

s'A/^([/)to+ — dm^{s'Un+) — dn^{s'Um+) on e. 
Furthermore, from (|6.6p we have 

dm+{s'Un+) = -m\n^be ^ 

on e, and, using s'Uf_m- = s'Uf^m+ = 0, we obtain 

/ s'Uf_m^ , s'Uf_n^ 

dn^{s Um+) = -be = -m^n^be . 

11^ 



It follows that 



s Af^{U)m+ = ■m\_n_be{ ^ j- ) on e. 

/i-i- h— 
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We have therefore shown that s'Af{U)m = on all edges if and only if, 

s'Uf ri_ s'Uf,n+ 

on all edges of K. 

Finally, we consider ni' Af{U)m = 2dm{'m'Un) ~ dn{m'Um). Using the fact that 
both m'_U f_m- and m'^U f^m+ vanish on e, we obtain from (|6.6p that, on e. 



h+ 

so 

, I n'_Uf_n^ m',Uf,n+ m'_Uf_n- 
m_)_A/^ ([/)m+ = m^n^be[m^n^ 2( m_^_m- )J. 

Hence, m'_^Af^(U)m+ vanishes on e if and only if. 

n_Uf_n- = — ■ ( ; m,m- ; ) on e. 

Note that this condition is not symmetric in /_ and /-|_. We thus obtain two 
conditions for each edge e. 

Since Uf G pI_^{K;S), with QnUfQn = on /, it follows that Uf is uniquely 
determined by the vector field Vf := UfU £ Vksif',^^)- The analysis above shows 
that U = J2f f^fUf ^ ■N'k.dK if and only if these vector fields satisfy 

(A) = t^^l± on e, and 

(B) n'_Vf_ = , ( 7 ^ - m'm^ — r^^) e, 

whenever an edge e is shared by faces /_ and /_(-. Therefore, there is an isomorphism 
between Mk,dK and 

(6.8) { {vf) e Yl n-3(/;M') I the Vf satisfy (A) and (B) }. 

To compute the dimension of the space (|6.8p we consider the relations (A) and (B) 
at a fixed vertex x of K. If (vf) is an element of the space (|6.8p define z S M'^ by 

(6.9) s'z = hfs'vf{x) 

for s chosen as tangents to each edge e meeting a;, and where / is a face meeting 
e. Note that the vector z is well-defined as a consequence of condition (A), and 
that for each face / containing x we have hft'vf[x) = t' z for all vectors t which are 
tangential to the face /. Using the expansion n_ = [771+ — (rn'j^m^ — m-\/m'_^n- 
we can then rewrite condition (B) at the vertex x as 

(6.10) hfn'j:Vf(x) = 2nfZ. 

From this discussion we can conclude that the dimension of the space (|6.8p . and 
hence dimAfk.dK ^ is at least 6(fc^ — 6fc + 10). To see this observe that 

dim[]Pfc_3(/,K') - 6(fc - 2)(fc - 1). 
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Furthermore, the conditions (A) and (B) represent a total of 6-3 • (A: — 4) = 18(fc — 4) 
constraints in the interior of the edges and 4 • 6 = 24 constraints at the vertices. 
Since 6{k - 2)(fc - 1) - 18(fc - 4) - 24 = 6{P - 6fc + 10), this is a lower bound for 
dim7\4,aif ■ 

We complete the proof by showing that elements of the space (|6.8[) are determined 
by 6(A:^ — 6fc + 10) degrees of freedom, in fact by degrees of freedom corresponding 
to the space 

To see this, for each vertex x pick a vector z ~ z{x) E M.^ and choose Vf{x) such 
that the relations (|6.9[) and (|6.10[) hold for all faces meeting x. This determines the 
vectors Vf{x) for all vertices x G df. We then define hfs'vf on each edge by the 
standard interior degrees of freedom, and m'vf with respect to both faces meeting 
e are determined similarly. These degrees of freedom on the edges correspond 
to the space J^^ 7-'fe_5(e; R^^). The normal components n'j^Vf are determined on 
each edge by condition (B). Finally, we must apply the interior degrees of freedom 
to Vf on /. We conclude that elements of the space Mk,dK are determined by 
12 + 18dimPfe_5(e) + 12 dim7'fc-6(/) = 6(fc2 - 6fc + 10) degrees of freedom. □ 

For U G Mk,dK and / G /S.2{K), A.f{U) is a polynomial vanishing on df, and so 
the quotient Af{U)/bf is a polynomial. We define Tf : Afk ok Vk-iif iQi3Qn) 

by 

Tf{U) = -hjK}{U)/hf. 

Lemma 6.4. // f-, f+ G A2{K), e = /_ n /+, and s is a unit vector parallel to e, 
then 

s'Tf^{U)s = s'Tf_{U)s on e, U £ Nk.dK- 
Proof. First we show that 

(6.11) dm+s'hf+{U)s^ drugs' Kf_{U)s one. 

Recall that dm+bf_^ = —m'_^n-be/h-, dm^bf_ = —m'_n+be/h+ on e and m'^n^ = 
m'_n+. We have on an edge e, dms'Af{U)s = 2dsdms'Un — dmdns'Us. Using (|6.6p 
we obtain 

dm+s'Un+ = -be—^^ s'Uf^n+, 

which is symmetric in /_ and /+ as a consequence of (j6.7|) and m'_^n^ = 77i'_77+. 

The identity (|6.1ip will follow if we show that dm+dn^s'Us — dm-dn_s'Us. 
Consider first the term 

V = bf_Uf_+bf^Uf^. 

Since QnUfQn = for / = /_, /+ and grad6/_ = — -^77,+ on /-). we derive that at 
the edge e 

1 TTl! JIj 
dm+dn+s'Vs = -be[ — sdm+Uf_S ^ 7-— S dn+U f_^_s\ 

/l-j- fl — 

= -{m'^n^)be[^s'dn_Uf_s + -^s'dn+Uf^s], 

and this expression is symmetric in /_ and /+ . Finally, consider terms of the form 
W — bfUf, where / is neither /_ nor /+. In this case 6/ = A/_A/^A, where A is 
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the barycentic coordinate associated the fourth face (7^ /, f-, f+) of K, and on e 
we have 

dm+dn+S Ws = As UfS. 

This is again symmetric in /_ and /+. Wc have therefore estabhshed (|6.1ip . 
Now, by definition, h+Af^{U) = -bf^Tf^{U). Therefore 

{m'^n^)-^h+h^d,n+s'Af^{U)s = bes'Tf^{U)s on e. 

By (|6.1ip . the left-hand side is unchanged if we interchange the subscripts + and 
— , so the same must be true of the right-hand side. □ 

Lemma 6.5. Let (Tf) G n/eA2(-R') Qn^Qn) be such that s'Tf S = s'Tf^s on 

e, whenever e = /_ n /+, /-, /+ € A2{K). Then there exist an S G Vk{K;S) such 
that QnSQn = Tf for all f e A2iK). 

Proof. We win define S S VkiK', S) by first specifying its vertex values, then speci- 
fying its moments of degree at most fc — 2 on the edges, then its moments of degree 
at most fc — 3 on faces, and then the moments of degree at most fc — 4 over the 
interior of K . 

Let a; be a vertex. We define the matrix S{x) S § by specifying the values 
s[S{x)sj where the Si are the tangents to the edges meeting at x (and so the Si 
form a basis for M.^). Namely we take s'^S{x)sj = s[TfSj with / the face containing 
Ci and ej. If i = j there are two possible choices of the face /, but they give the 
same result by assumption. 

For the interior degrees of freedom on an edge e we use the basis s, m_, m+ of 
R^, and let e 7'fc(e;S) be given by 

S leS = S 1 f_S = S 1 fj_S, s lem- = s 1 f_m-, s lem+ = s 1 f^m+, 
Then we define 5|e by 



f{S^T,)Vds = 0, V&Vk-2{e;l 

J e 



Similarly, for the interior degrees of freedom on each face we let QnSQn inherit the 
moments from Tf, while the data for SPn is taken to be zero. The interior degrees 
of freedom on K are all taken to be zero. □ 

As a consequence of the two previous lemmas, there is a map 

(6.12) A4,aA' ^ Vk^i{K- §), U ^ 5(C/), 

such that QnSiU)Qn ^Tf{U) for all faces f £ A2{K). 

We are finally ready to compute the dimension of the space Afk defined in (|6.ip . 

Theorem 6.6. For fc > 3, the dimension of the space Afk is k{k'^ — 6fc + 11). 

Proof. Let S G A/fe- By Lemma [6.21 S must be zero on each edge and so can be 
written S ^J2f bfSf + bKSK, where Sf G Vl_^{K] S) and Sk G Vk-i{K; §). Now 
ef{bKSKnf) vanishes on / since bK does, while 

bf 

dnibKQnSxQn) = {dnbK)QnSKQn = --;—QnSKQn OU /. 
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Thus 

(6.13) AfibKSK) = ^QuSkQu- 

In particular, Kf{hKSK) vanishes on df . It foUows that if S* e Nk and we define 
U = bfSf then U S Numk- Therefore, the map ([/, Sk) ^ U + IkSk defines 
an isomorphism from 

{ (f/, Sk) e Mk,aK x Vk-A{K; §) | QuSkQu = Tf{U) on each face / } 

onto Nk ■ 

Finally, note that a matrix field of the form bxY, V G Vk-iiK-jS), belongs to 
A/fc if and only if F G N^^^. Therefore, using the map (|6.12p . the mapping 

^k,aK X ^ Nk, {U, V) ^ {U, S{U) + V), 

is an isomorphism. It follows that dim Nk = dim Nk.dK + dim N^^^ and using 
Lemma 15751 and Lemma lOl wc get dim7\4 = 6(fc^ — 6fc + 10) + (fc — 3)(fc — 4)(fc — 5) = 
fc(fc2-6fc + ll). □ 

7. The space oe divergence-free matrix fields with vanishing normal 

TRACES 

Recall that the space 

Mk=^MkiK)^{SerkiK;S)\divS^O onK, P„5 = 0, /eAalA')} 

appears in the degrees of freedom for the finite element space S^, C i7(div, fi; §) 
introduced in Sections[3]and[4l Therefore, a derivation of the dimension of this space 
is fundamental for our theory, while a construction of a (dual) basis for the space 
Mk is necessary for the implementation of the method. The dimension formula 
will be a simple consequence of the following lemma, in which 'P^^^{K; R'^) := {v G 
Vk+3{K;R^)\v = ondK} = bKVk-i{K;B?). 

Lemma 7.1. (1) The operator curl curl* maps Nk+2{K) onto Nlk{K). 

(2) {TeNk+2\ curlcurl*r = 0} = e[pO+3(A';R3)]. 

(3) The following sequence is exact: 

(7.1) ^ P°+3iK;R')^Nk+2{K) Mk{K) ^ 0. 

Proof. It follows directly from (|5.7p and (|5.14p that curl curl* 7\4+2 C Mk- Hence, 
to prove the first statement we need only show that A4k C curl curl* Nk+2 ■ Let 
S € Mk- Since div5 = 0, it follows from the exactness of the complex (|2.6p that 
there is a T G Vk+2{K\ §) such that S ~ curl curl* T. The proof will be completed 
by constructing a vector field u G Vk+7i{K\lS.^) such that 

(7.2) g„(T-e(w))g„ =0, A/(r-e(it)) =0, on each face /. 
Note that since S & Nik it follows from (|5.7p that 

rot/ rot^ QnTQn = P„ (curl curl* T)P„ = PnSPn = on each face /. 

Hence, from the exact sequence (|5.9p we conclude that for each face / G A2{K) 
there is a vector field Vf G Vk+sif',^^), with PnVf = 0, such that QnTQn = 
tf{vf). The vector fields Vf are uniquely determined up to a 2D rigid motion 
and hence we may normalize them so that / s'vfds = on each edge e G f. 
Since ds{s'vyj_) = s'Ts = ds{s'vf^) on each edge, we obtain that PsVf_ = PgVf^ 
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on each edge e = /+ n /-. As a consequence, there is u e Vk+siK-jM."^) such 
that QnvQn = Vf on each face /. Then Q„e(w)Qii = = QnTQn, i.e., 

QnUQn = on each face /, where U — T — e{v). This implies, in particular, that 
Us and grad(s't/s) vanish on each edge e £ Ai{K). Therefore, 

(7.3) s'Af{U)s = 2ds{n'Us) - dn{s'Us) = on 9/ 
for each face /. 

Next, observe that, by ([5l^ . rot/ A/(C/) = CnSPn = 0. Hence, ([ETU]) implies 
that there is a scalar field qf £ Vu+iif]^), uniquely determined up to a linear 
function on /, such that Af{U) = grady gradj qf. On each edge e S 9/, we have by 
(|7.3p that = s'Af{U)s = 9,(7/. It follows that we can assume that = on 9/. 
Hence, there exists another vector field w in Pk+3{K;M.^) such that = and 
n'w = qf on each face. Recall by (|5.13[) that A/(e(w)) = grady gradj qf ~ Af{U). 
Hence, if we let u = w + w, then the relation (|7.2|) holds. This proves the first 
statement. 

We now prove the second statement, li T = e{v) for some v G V^_^^{K;M.'^), 
then curl curl* T = 0, and, by ([5l^ . 

(7.4) QnTQn ^efiQnv), Af{T)^grSidfgY&d*f{n'v) on/, 

for each face /. Since v vanishes on /, the right hand sides of these equations vanish, 
and so T belongs to Mk+2{K). Conversely, if T S Nk+2{K), then, by the exactness 
of the sequence (j2.6|) . T ~ e{v) for some v G Vk+3iK;M.^), which is determined 
uniquely if we require that 

(7.5) / s'vds = 0, e G Ai(/v). 

J e 

(The functionals v i-^ J^s'vds, e £ Ai{K), form a set of degrees of freedom for T, 
the null space of e.) From the first equation in (|7.4p and (|7.5p . we find that QnV 
vanishes on each face /. Therefore the entire vector v vanishes on each edge e. 
Using the second equation in (|7.4[1 . we see that n'v vanishes on each face as well, 
so V e V^^^{K;M.'^). This completes the proof of the second statement. 

The third statement is an immediate consequence of the first two and the fact 
that Tn7'^+3(i^;M3) = 0. □ 

Theorem 7.2. Fork > 4 the space Mk{K) has dimension (fc + 2)(A: - 2)(fc - 3)/2. 

Proof. Using first the short exact sequence in the lemma and then the dimension 
formula in Theorem 16.61 we get 

dimA^fc = dim7Vfc+2 - dime[7'j?+3(A';]R^)] 

= {k + 2){k^ - 2fc + 3) - dimPk-iiK;M.^) = {k + 2){k - 2)(A:-3)/2. □ 

To conclude this section, we construct a basis for the space AAk{K) for fc = 4 and 
= 5. (Alternatively a basis could be constructed for any k using computational 
algebra software.) For this we use the following lemma, similar to Lemma |7. II 

Lemma 7.3. (1) The operator curl curl* maps bKAfj^ ^i^"^) onto A4k{K)- 

(2) { T e feA-A/'^2 I curl curl* T = } = e[blPk-5iK; M')]. 

(3) The following sequence is exact: 

(7.6) ^ blPkMK;R') ^ bKM^_,{K) MkiK) ^ 0. 
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Proof. Note that 6xAA°_2 ^ by ([6l3l) . and so curl curl* 6/fA/'"_2 

First we prove 2. Let w £ 7-'fc_5(/C; M"^). By the Leibniz rule 

e{bjiw) = bj^ e{w) + bK[{gr&dbK)w' + wigr&dbK)']- 
Clearly bx £{w) G ^"^d, recalling that gradfe^s- — —bfUf/hf, we see that 

(grad6if )u;' + zi;(grad6/f)' £ Thus €{bj^w) G bKM^_2^ giving the inclusion D. 

Conversely, if T e bkN^_2 with curl curl* T = 0, then, by Lemma [TTTl T = e{bKv) 
for some v G ; R"^) and we need to show that w = on dK. Using the 

Leibniz rule and the fact that T vanishes on dK, we get that nv' + vn' vanishes on 
each face. We conclude that v vanishes on the face, using the elementary identity 
w = (/ + Qn){nv' + vn')n/2. 
It follows that 

dim[curlcurl*(6/<-A/';;?_2)] = dimA/;"_2 ' dimPk-siK]^^) = dimTW^, 
where we have used Lemma 16.11 and Theorem 17.21 The exactness of (|7.6p , and so 
also the first statement of the lemma, follows. □ 

Thus for k = 4, curl curl* is injective on fe^AA^, and so a basis for A44 = 
curl curl* (fei^-A/J) is computable directly from a basis for A/J , which may be ob- 
tained directly from (|6.3[) . 

Now let fc = 5. The map curl curl* is not injective on bxAf^, but has a kernel 
of dimension 3. In this case, the representation (j6.3p presents an arbitrary element 
S e as 

eeAi{K) feA2{K) 

Fix a particular face /o G A2{K) and define TVj" as the subspace of 5 e A/J for 
which Sf„ = in this representation, clearly a subspace of codimension 3. We 
claim that curl curl* is injective on the space bxMi^, and hence a basis for A/fs 
can be computed from a corresponding basis of A/J". The injectivity follows since 
if w G M^, with e{b\w) £ bxM^^ , then we get w = arguing as in the proof of 
Lemma 17.31 

8. A Discrete Elasticity Complex 

The results of Section [7] above completes the description of the finite element 
spaces Vh and E/i, introduced in Sections [3] and IH and therefore also of the finite 
element method p. II) . However, as already indicated in the beginning of Section [SI 
there are more structures hidden in the construction above. In fact, the spaces Vh 
and Eft, are constituents of a discrete elasticity complex of the form 

(8.1) T^Wh^Qh ^^^^^^ E, ^Vh^O, 

where Wh C H^{il;M.'^) and 8^, C iJ(curlcurr, O; §) are piecewise polynomial 
spaces with respect to the triangulation Th. Furthermore, there exist interpolation 
operators n)^, 11^, 11^ and II]^ such that the diagram 

T C°°(r2;R3) ^ C°°(f7;§) > C°°{n;S) C°°{n;R^) > 

TP Tir <^ curl curl* „ div ^ 

T ^ Wh > Qh > E/i y Vh > 
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commutes. The spaces and Vh, and the associated interpolation operators 11^ 
and n]^, have been introduced above, so it remains to define the spaces and 
Qh, and the associated interpolation operators. 

The discrete complex (|8.ip can be defined for all polynomial levels, i.e., the two 
final spaces S/j and Vh can be taken as any of the pairs in the family introduced in 
Section m However, in order to simplify the discussion below we will only discuss 
the lowest order case introduced in Section [3l i.e., C H{div,rt;§) consists of 
piecewise quartic matrix fields with linear divergence, while Vh G is 
composed of piecewise linears. 

We will first describe the corresponding space Qh C i?(curlcurr, fi; §). Locally 
on each tetrahedron this space consists of functions in VeiK; §), which is a space of 
dimension 504. In order to specify the degrees of freedom on a tetrahedron K G Th 
we introduce the polynomial space 

Eif) = {ve Prif; K) \v\of^ 0, drnv ds^O^ee Ai(/) } 

for each face / G A2{K). The dimension of E{f) is 12. Also for each e G Ai{K) de- 
fine the operator : C°°{K;S) C°°(e;e-L) byre(5') = 2dsQsiSs) -grad^^ s'Ss. 
Note that if / G A2{K), with n = 7i/, and e G Ai(/), then nTe(S') = s'Af{S)s on 
e. 

The 504 degrees of freedom used to define the finite element space Qh are the 
following: 

(1) S and curl curl* S at each vertex, 48 degrees of freedom, 

(2) s'Ss vds,v e Viie; R), e G Ai(/-i:), 30 degrees of freedom, 

(3) jlTe{S) -vds, V G r5{e;QsM.^), e G Ai{K), 72 degrees of freedom, 

(4) (curl curl* S*) : W ds, We 7^2(6; Q^S), e G Ai{K), 90 degrees of 
freedom, 

(5) rotgj_ PgSQs ds, e G Ai{K), 6 degrees of freedom, 

(6) ])QnSQn ■■ ef{v)dxf, V G &/7'4(/;Q«M^), / G A2{K), 120 degrees of 
freedom, 

(7) J^is'Af{S)m)vds, V G Vi{e;R) with v ds = 0, / G A2(K), e G Ai(/), 12 
degrees of freedom, 

(8) J^Af{S)dxf, / G A2{K), 12 degrees of freedom, 

(9) J^.Af{S) : gradj grad} v dxj, v G E{f), f G A2{K), 48 degrees of freedom, 

(10) curl curl* S -.Tdx^T e Mi{K), 6 degrees of freedom. 

(11) 5* : e{v) dx, v G r^{K; R^), 60 degrees of freedom. 

These degrees of freedom are unisolvent for the space Ve{K;§). This is in fact a 
consequence of the following result. 

Lemma 8.1. Let f G A2(/^) be fixed and assume that S G Ve{K;§) with all the 
degrees of freedom (l)-(9) associated to all g G A(/) equal to zero. Then Qn SQn 
and Af{S) are identically zero on the face f. 

Proof. Let / G A2(A') be fixed, and asssumc that S G Ve{K;§) has all degrees of 
freedom associated with the subsimplexcs of / equal to zero. We start by observing 
that the degrees of freedom (l)-(4) implies that 

(8.3) s'Ss = 0, curl curl* 5 = 0, re(S') = 0, on e G Ai(/). 
Next we will show that 



(8.4) 



rot/ rot J- QnSQn ~ P„ (curl curl* 5)P„ = 
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on /. Note that by ()8.3p this quantity vanishes on df, i.e., rotf rotj QnSQn G 
bfri{f;M.Pn). Furthermore, using ([5l^ . for each U £ 'Pi(/;RP„) 

/ rotfiot*fQnSQn ■.Udxf = - / [rot/ (5s) • (Un) + {Ss) ■ curl f{Un)] ds. 
Jf Jdf 

Also, by (|5.ip we obtain rotj(S's) = — [5s(s'C„S's) + 9m(m'C„S's)]n. However, 
the fact that m'TeiS) = on e G Ai(/) imphes that dms'Ss = 2dsm'Ss or 
dmm'CnSs = —2dss'CnSs. Therefore, we can conclude that rot/ (5*8) = ds{s'CnSs) 
on e G Ai(/). On the other hand, since s'Ss = on e we have from (|5.3p that 

(S's) • curl/([/?i) = -(S-s) • [C„grad(n'[/7i)] 

= (CnSs) ■ grad(n'[/n) = (s'C„S's)(9,(n';77i). 

Since S is zero at each vertex of / we can therefore conclude that 

(8.5) ( T0ifmt}QnSQn-U dxf = - j ds[[s' C„Ss){n'Un)]ds = Q 

J f Jdf 

for all U G 7'i(/;RP„), and hence (|8.4p follows. As a consequence, we obtain from 
()5.9p that QnSQn = ^f{u) for a suitable u G 'P7(/; (5„M^), where u can be chosen 
such that s'uds = for each edge e G Ai(/). Therefore, since s'S's = ds{s'u) = 
on each e G Ai(/), we conclude that s'w = on 9/, and that u = at each vertex. 
Furthermore, on each edge e G Ai(/) 

= m'TeiS) = 2ds{m'Ss) - d,n{s'Ss) = dl{m'u), 

and therefore w = on df . Hence, u vanishes on / by the degrees of freedom (6), 
and so does QnSQn- 

Next we will show that rot f Af{S) = on /. By (|5.14p we have that rot /A/ (5) ~ 
C„(curlcurl* S)Pn, and therefore (|8.3[) implies that rot/ A/(5) is zero on df. Hence, 
it is enough to show that 

(8.6) [ TotfAfiS) -.Vdxf = 0, l^GPi(/;g„§P„). 
Furthermore, by the degrees of freedom (8) we obtain that 

[ TotfAfiS) -.Vdxf = [ [Af{S)s]-iVn)ds^ [ [m'Af{S)s]-{m'Vn)ds, 
Jf Jdf Jdf 

for any V G Vi{f] QnSPn), where the last identity holds since s'Af{S)s = n'TeiS) = 
on e G Ai(/). From (7) above, this will be zero if wc can show that 7Ji'Af{S)s ds 
vanishes for all e G Ai(/). However, ■m'Af{S)s ~ ds{m' Sn) + dm{s' Sn) — dn{s' Sm) . 
Therefore the degrees of freedom (5) imply that J 'm'Af{S)sds = J ds{m'Sn)ds = 
0. We therefore conclude that rot/ A/ (5) is zero on /, and from (|5.10p wc can 
conclude that A/(S') = grad/grad/W for a suitable scalar field v G Vjif^M.), 
where we assume that v is chosen to be zero at each vertex. Furthermore, the 
fact that = ri're(5) — s'Af{S)s = d^v implies that w = on df. In par- 
ticular. grad/U = at each vertex. Therefore, v G E{f) since (7) implies that 
— s{s' Af{S)m) ds = sdsdmV ds ~ ~ dmV ds, and, as a consequence of (9), 
V = 0, and so is A/(5'). □ 

The lemma above implies that the set of functionals (l)-(ll) are unisolvent for 
the space VeiKiS). This follows since if all the functionals (l)-(ll) are zero then 
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S g A/e by this lemma, and hence cmdcmd* S G Ma- By the degrees of freedom 
(10) cmdcmd* S ~ 0, and by Lemma FfTTl and the degrees of freedom (11), 5 = 0. 

The finite element space Qh is defined as all functions which belongs to Pe{K; S) 
for all K £ Th and with the continuity conditions induced by the degrees of free- 
dom. Hence, as a consequence of Lemma |8.1[ the variables QnSQn and A/(S') are 
continuous for all 5 G 8/i and / G A2(7/i), and therefore Theorem 15.21 implies that 
Qh C iJ(curlcurr, ri; §). The degrees of freedom (l)-(ll) also defines a canonical 
interpolation operator : C°°(f2; S) ^ Qh by requiring that it reproduces all the 
functionals. 

Finally, we need to describe the finite element space Wh C H^{Q.]M.^). Locally 
on each K G 7/i this space is taken to be VjiK; R'^), which is a space of dimension 
360. A vector field w G 'P7{K;M.'^) is uniquely determined by: 

(1) the values of w and its first order derivatives at each vertex, 48 degrees of 
freedom, 

(2) J^w qds, q <E 'P3{e;R^), e G Ai(A'), 72 degrees of freedom, 

(3) J^efiQnw) : ef{q)dxf, q G bfP^iflQnM.''), f G A2{K), 120 degrees of 
freedom, 

(4) drn{n'w) ds, f G A2(A'), e G Ai(/), 12 degrees of freedom 

(5) grady grad}(n'w) : grady grad} (g) dxj, q G E{f), f G A2(A:), 48 degrees 
of freedom, 

(6) Jj^ e{w) ■ e{q) dx for all q G Vj{K; IS?), 60 degrees of freedom. 

It is rather straightforward to check that this set of functionals is unisolvent for the 
space V7{K; M^). The finite element space Wh is then defined as the set of functions 
which are locally in Vt{K; M^) and with the continuity induced by these degrees of 
freedom. Hence, all w G Wh are continuous across faces / in A2(7;i), and they are 
C"^ at the vertices. The associated interpolation operator H)^ : C°°(17;M^) — > Wh 
will reproduce all degrees of freedom given by (l)-(6). 

By using the degrees of freedom for the discrete spaces, one may also verify that 
the diagram (j8.2p commutes. This verification is tedious, and we will drop most 
of the details. However, we will illustrate what needs to be done. For example, as 
part of the verification of the relation H|' o curl curl* = curl curl* oH^ we must show 
that for T ={I- n^)5, S G C°°(17; §), we have (cf. degrees of freedom (3) for Y^h) 

^ P„ (curl curl* r)P„ -.Udxf = j rot/ rot} QnTQn ■ Udxj = 0, 
J C„(curlcurl*r)P„ : Vdxj ^ j lotf kf{T) : Vdxj =0, 

for U G Pi(/,MP„), V G Pi(/,(3„§P„), and aU / G A2{Th). However, the ver- 
ifications of these identities are almost identical to arguments leading to the for- 
mulas (|8.5p and (|8.6p in the proof of Lemma 18.11 above. The two identities above 
show that the degrees of freedom (3) for E/j are zero for functions of the form 
(H^ curl curl* — curl curl* H^)^ G S;, . Similar, but simpler, arguments can be used 
to show that the other degrees of freedom for S/j are zero. Hence, the identity 
H^ o curl curl* = curl curl* oH® follows, and corresponding verifications can be 
done for the identity H® o e = e oH)^. 
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